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Abstract —The large-scale structure of the universe is comprised of virialized blob-like clusters, linear filaments, sheet-like walls and 
huge near empty three-dimensional voids. Characterizing the large scale universe is essential to our understanding of the formation 
and evolution of galaxies. The density range of clusters, walls and voids are relatively well separated, when compared to filaments, 
which span a relatively larger range. The large scale filamentary network thus forms an intricate part of the cosmic web. 

In this paper, we describe Felix, a topology based framework for visual exploration of filaments in the cosmic web. The filamentary 
structure is represented by the ascending manifold geometry of the 2-saddles in the Morse-Smale complex of the density field. We 
generate a hierarchy of Morse-Smale complexes and query for filaments based on the density ranges at the end points of the filaments. 
The query is processed efficiently over the entire hierarchical Morse-Smale complex, allowing for interactive visualization. 

We apply Felix to computer simulations based on the heuristic Voronoi kinematic model and the standard ACDM cosmology, and 
demonstrate its usefulness through two case studies. First, we extract cosmic filaments within and across cluster like regions in Voronoi 
kinematic simulation datasets. We demonstrate that we produce similar results to existing structure finders. Second, we extract different 
classes of filaments based on their density characteristics from the ACDM simulation datasets. Filaments that form the spine of the 
cosmic web, which exist in high density regions in the current epoch, are isolated using Felix. Also, filaments present in void-like regions 
are isolated and visualized. These filamentary structures are often over shadowed by higher density range filaments and are not easily 
characterizable and extractable using other filament extraction methodologies. 

Index Terms —Morse-Smale complexes, tessellations, cosmology theory, cosmic web, large-scale structure of the universe. 
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1 Introduction 

t scales from a megaparsec to a few hundred mega¬ 
parsecs 1 , the universe has a web-like appearance. In 
the cosmic web [1], [2], galaxies, intergalactic gas, and dark 
matter have aggregated in an intricate wispy spatial pat¬ 
tern marked by dense compact clusters , elongated filaments 
and sheetlike walls , and large near-empty void regions. The 
filaments, stretching out as giant tentacles from the dense 
cluster nodes, serve as transport channels along which mass 
flows towards the clusters. They surround the flattened walls, 
which are tenuous, membrane-like features in the cosmic mass 
distribution. 

All structures and objects in the universe emerged out 
of primordial fluctuations that were generated during the 
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1. A parsec is the standard unit of measurement of distances in the cosmos. 
A parsec is 3.26 times the light-year, the distance light covers in a year. A 
megaparsec is a million parsecs, the typical scale of measurement of size of 
the large scale structures in the universe. 


inflationary era , moments after its birth, as the universe 
underwent a rapid phase of expansion [3], [4], [5]. The 
quantum fluctuations generated during this phase manifest 
themselves as fluctuations in the temperature of the cosmic 
microwave background [6], [7], [8]. The gravitational growth 
of these density and velocity perturbations has resulted in the 
wealth of structure that we see in the Universe. The web¬ 
like patterns mark the transition phase from the primordial 
Gaussian random field to highly nonlinear structures that have 
fully collapsed into halos and galaxies. As our insight into the 
complex structural pattern of the cosmic web has increased 
rapidly over the past years, it has become clear that the 
cosmic-web contains a wealth of information on a range of 
cosmological and astronomical aspects and processes. 

An important illustration of the cosmological significance of 
the cosmic web concerns its dependence on the nature of dark 
energy and matter, the dominant but as yet unidentified forms 
of energy and matter in the Universe. One telling example of 
this is the recent realization that cosmic voids are sensitive and 
useful probes of the nature of dark energy and dark matter and 
testing grounds for modified gravity theories [9], [10], [11], 
[12], [13], [14], [15], [16]. As the cosmic web is first and 
foremost defined and shaped by the gravitationally dominant 
dark matter, it would be of considerable importance to be 
able to obtain detailed maps of dark matter distribution. In 
recent years, great strides have been made towards this goal 
as gravitational lensing of distant galaxies and objects by the 
dark matter have enabled an increasingly accurate view of 
its spatial distribution [17], [18]. Initial efforts concentrated 
on the detection and mapping of the deep potential wells 




IEEE TRANSACTIONS ON VISUALIZATION AND COMPUTER GRAPHICS, VOL. X, NO. X, JANUARY XXXX 


2 


of the nodes in the cosmic web, i.e., of galaxy clusters. 
Recent results have opened the path towards the mapping 
of filaments via their lensing effect on background sources 
[19]. The identification of the structural components of the 
cosmic web is also important for our understanding of the 
relation between the formation, evolution, and properties of 
galaxies and the structural environment of the cosmic web. A 
direct manifestation of this is the generation of the angular 
momentum of galaxies. This is a product of the torqueing by 
the large-scale tidal force field [20], [21], [22]. While these 
force fields are also the agents for the formation and shaping 
of filaments, we would expect that this results in the alignment 
of the spin axis of galaxies with respect to cosmic filaments 
[23], [24], [25]. 

The identification, description, and characterization of the 
elements of the cosmic web is a non-trivial problem. Several 
characteristics of the mass distribution in the cosmic web have 
made it an extremely challenging task to devise an appropriate 
recipe for identifying them: 

a) The cosmic web is a complex spatial pattern of connected 
structures displaying a rich geometry with multiple morpholo¬ 
gies and shapes. 

b) There are no well-defined structural objects at a single 
spatial scale or within a specific density range. Instead, ele¬ 
ments of cosmic web are found at a wide range of densities 
and spatial scales. This is a consequence of the hierarchical 
evolution of structure formation in the universe, such that 
smaller high-density structures merge to form larger objects. 

c) There is a clear anisotropy in the structures of the 
cosmic web, a consequence of gravitational instability. The 
structures in the cosmic web exhibit elongated and flattened 
characteristics. 

The attempts to analyze the structure of the cosmic web 
have a long history. The absence of an objective and quan¬ 
titatively accurate procedure for identifying and isolating the 
components of the cosmic web has been a major obstacle in 
describing it. In recent years, more elaborate and advanced 
techniques have been developed to analyze and describe the 
structural patterns in the cosmic web. Nonetheless, a consensus 
on the proper definition of filaments is yet to be achieved. In 
the subsequent subsection 1.1 we present a short account of 
the available techniques and the definitions on which they are 
based. 

1.1 Related Work 

Statistical measures such as the auto-correlation function [26] 
of the matter distribution in the web have been the mainstay of 
cosmological studies over many decades. However, while this 
second-order measure of clustering does not contain any phase 
information (one may e.g. always reproduce a distribution with 
the same 2nd order moments and random Fourier phases), the 
auto-correlation function is not sensitive to the existence of 
complex spatial patterns. Higher order correlation functions 
only contain a very limited amount of such structural informa¬ 
tion, while in practical observational circumstances, it quickly 
becomes cumbersome to measure them as the magnitude of 
the error increases drastically with increasing order. 


The first attempts towards characterizing complex geometric 
patterns in the galaxy distribution mainly involved heuristic 
measures. Early examples of techniques addressing the global 
connectivity of structure in the Universe are percolation anal¬ 
ysis [27] and the minimum spanning tree of the spatial galaxy 
distribution [28], [29]. While these are useful global descrip¬ 
tions, they do not capture and describe local characteristics of 
the mass distribution. 

More elaborate and advanced techniques have been devel¬ 
oped in recent years. Several of these methods apply sophis¬ 
ticated mathematical and visualization techniques, involving 
geometric and topological properties of the cosmic mass 
distribution. There are a multitude of different methods for 
detecting filaments, based on a range of different techniques. 
We may recognize several categories of techniques. 

One class of methods seeks to describe the local geometry 
on the basis of the Hessian of the density field [30], [31], [32], 

[33] or closely related quantities such as the tidal force field 

[34] , [35] or the velocity shear field [36], [33]. The Hessian 
provides direct information on the local shape and dynamical 
impact of the corresponding field. The morphological elements 
of the cosmic web are identified by connecting the areas within 
which a specific range of anisotropies is registered. 

These studies concentrate on a single scale by appropri¬ 
ately smoothing the field, and do not consider the multi¬ 
scale nature of the cosmic mass distribution. The Hessian 
based Nexus/MMF technique [31], which was perfected into 
a versatile and parameter-free method [33], implicitly takes 
into account the multi-scale nature of the web-related fields. 
It accomplishes this by a scale-space analysis of the fields. 
At each location the optimal morphological signal is extracted 
via the application of a sophisticated filter bank applied to 
the Hessian of the corresponding fields in scale space. The 
application of this machinery has enabled thorough studies of 
the hierarchical evolution and buildup of the cosmic web [37]. 

A promising and highly interesting recent development has 
opened up the path towards dynamical analysis of the evolving 
mass distribution in full six-dimensional phase-space (in which 
the position of each mass element is specified by its space 
coordinates and velocity/momentum). In the 6D phase space, 
the cosmic mass distributions defines a 3D sheet. Indepen¬ 
dently, three groups arrived at tessellation based formalisms 
that exploit the evolving structure and folding of the phase 
space sheet in phase space [38], [39], [40] (also see e.g. 
[41]). The number of folds of the phase space sheet at a 
given location indicates the number of local velocity streams, 
and forms a direct indication of the morphology of the local 
structure. Interestingly, the resulting characterization of the 
web-like distribution, the Origami formalism of Neyrinck [40] 
for example, appears to resemble that of the Nexus/MMF 
formalism [37]. 

An entirely different class of techniques, based on statistical 
methods, have also been used to recover the filamentary 
patterns in the Universe. The key idea behind these techniques 
is to treat the galaxy distribution as a Markov point process. 
Within this class, particularly worth mentioning is the Bisous 
or the Candy model, which has been used by Stoica et al. [42] 
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to detect the filamentary network 2 . In this model, one places 
a random configuration of interacting geometric cylindrical 
objects on the point process to detect filamentary structure. 
It has been developed into a versatile, statistically solid yet 
computationally challenging formalism for the identification 
of filaments in a spatial point distribution, such as N-body 
simulations and galaxy redshift surveys [44]. An additional 
example of a method involving statistical analysis of a ge¬ 
ometric model is that of Genovese et al [45], which seeks 
to describe the filamentary patterns of the cosmic web in a 
non-parametric way by recovering the medial axis [46] of the 
point-set of galaxies. 

The fourth major class of methods, the one which we will 
also pursue in this paper, exploits the topological structure of 
the cosmic mass distribution. The fundamental basis of these 
methods is Morse theory [47]. The geometric structure of the 
Morse-Smale complex [48] naturally delineates the various 
morphological components on the basis of the connections 
between the critical points of the density fields and the higher 
dimensional cells that are incident on the critical points. 
Various Morse theory based formalisms have been applied 
to the identification of components of the cosmic web. One 
of the first applications concerned the detection of voids in 
the cosmic density field. The Watershed Void Finder [49] 
identifies these with the watershed basins around the density 
minima. The Spine Web procedure [50] extended the watershed 
transform towards the detection of the full array of structural 
components, filaments, walls and voids. These techniques use 
a user-defined filter to incorporate the multi-scale structure of 
the cosmic density field. 

A natural topological means to address the multi-scale 
topological structure emanates from the concept of persistence 
[51]. It provides a natural recipe for detecting and quantifying 
the components of the cosmic web in a truly hierarchical 
fashion. Sousbie [52], [53] has exploited and framed this in an 
elegant and impressive framework, the DisPerSE formalism. 
Following the construction of the Morse-Smale complex, 
they proceed to simplify it. The simplification proceeds by 
canceling pairs of critical points iteratively, where each pair 
represents a structure in the cosmic web. Topological persis¬ 
tence is invoked to order the critical point pairs. However, this 
measure of importance is not unique, and one may consider 
alternatives, dependent on the specific interest and purpose. 

In effect, to tackle similar issues in other visualization areas, 
a range of variations have been proposed in other studies 
[54], [55], [56]. Weinkauf et al [56] describe the concept of 
separatrix persistence, where they compute the strength of sep¬ 
aration of points on a separatrix curve (in 2D) connected to a 
saddle as the sum of the absolute differences of function values 
of the saddle and the extrema connected to it. This concept 
is extended to 3D separating sheets by Gunther et al. [54]. 
Reininghaus et al [55] develop the concept of scale-space 
persistence where they accumulate the absolute difference in 
function value measure of critical points across a hierarchy of 
derived functions. The set of derived functions are generated 

2. See Stoica et al. [43], for a detailed description of this object point 
process which is used to characterize an observed point distribution by a 
fitting procedure using global optimization techniques. 


by smoothing the function using a family of Gaussian kernels 
of increasing variances. This is similar to the Multi-scale 
Morphology Filter Nexus/MMF [31], [33] described above. 
Both methods adopt the scale space formalism as the first 
step to detect features at multiple scales. However, scale 
space persistence and separatrix persistence, disregard specific 
density regimes of interest and are potentially inappropriate 
when small scale features with specific density characteristics 
are of interest. 

1.2 Present study: contributions 

In the present study, we describe and introduce a technique 
for the identification of filaments based on the topological 
characteristics of the density field. A key aspect of the pro¬ 
posed technique is its interactive nature, involving a tunable 
density parameter. Specifically, we describe the following 
contributions: 

a) We describe Felix 3 : a topology based framework for 
visual exploration of filaments in the cosmic web. In particular, 
we develop a query framework to extract filamentary structures 
from a hierarchy of Morse-Smale complexes of the density 
field. The filaments in Felix are parameterized by the density 
values of the maxima and the 2-saddles that define them. 

b) Using Felix, we develop a semi-automatic structure 
finder that classifies galaxies as cluster/filamentary or not. 
We demonstrate its efficiency through two tests. First, using 
the Voronoi Kinematic model as a benchmark, we demon¬ 
strate that we are able to recover the classification with 
high efficiency. Second, we show that the classifications are 
quantitatively comparable to, and in several cases better than, 
existing classifiers. 

c) We investigate the nature of filaments in two different 
density regimes from the ACDM simulations. The first con¬ 
cerns filaments in the high density regions around compact 
dense clusters, which are known to function as the transport 
channels along which matter moves into the clusters. A second 
regime concerns the tenuous low-density filaments found in 
low-density void regions. In the supplemental document, we 
present an additional experiment, where we investigate the 
nature of three classes of filaments in a relatively cleaner 
region of a ACDM dataset. 

d) In the supplemental document, we describe an efficient 
structure based volume rendering enhancement routine that 
allows us to highlight the density distribution in regions that 
are close to the selected features. 

The distinction between noise and significant structures is 
often ill-defined, and at occasions noise may be confused 
with genuine structures in the hierarchically evolved mass 
distribution (see Figure 4 and the caption thereof for an illus¬ 
tration). This problem is more pronounced when one studies 
the properties of tenuous filaments and walls in low density 
void-like regions. For the understanding of the formation and 
evolution of galaxies in such regions, we need to assess 
the possible dependence of galaxy and halo properties on 
the morphology and density of the local environment. This 
must be based on the successful extraction of filaments in 

3. The name Felix is formed from an abbreviation of /T/ament explorer. 
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Fig. 1: Ascending manifolds of a 2-saddle (yellow sphere). 
The scalar function is a sum of two 3D Gaussians centered 
on either side of the volume. The two arcs incident on the 2- 
saddle constitute the ascending manifold and terminate at the 
two maxima (red spheres) of the scalar function. 

low density regions and the correct identification of galaxies 
associated with them. In view of this, we include an interactive 
handle on the density regimes so that one can concentrate on 
and probe structures in specific density regimes. 

The remainder of the paper is organized as follows. Sec¬ 
tion 2 introduces necessary background material. Section 3 
describes Felix. Section 4 introduces the cosmological datasets 
used in the experiments. Section 5 discusses the application of 
Felix as a structure finder and as a tool for exploring filaments 
in different density regimes. Section 6 concludes the paper by 
summarizing the main results and possible future directions. 

2 Background 

This section reviews relevant background on Morse functions, 
the Morse-Smale complex, and topological simplification. This 
is a necessary prerequisite for understanding the definition of 
filaments and extraction methods described in the subsequent 
sections. 

2.1 Morse theory and the Morse-Smale complex 

Let / : M —» R be a real-valued scalar function defined on a 
manifold M. Critical points of / are points of / where the 
gradient of / vanishes i.e., V/ = 0. Morse theory is the study 
of the relationship between the topology of level sets of scalar 
functions and the critical points of the function. The function 
/ is said to be a Morse function if all of its critical points 
are non-degenerate i.e., the Hessian of /, equal to the matrix 
of second order partial derivatives, is non-singular. The non¬ 
degeneracy condition imposes a locally quadratic form for / 
within a small neighborhood of its critical points. In other 
words, using a coordinate transformation, the function near 
a critical point p of the ^-dimensional manifold M can be 
written as f p (x) = f(p) =b x\ ±x% ± ... ± x*. The index of p is 
equal to the number of negative quadratic terms in the above 
expression. In 3D, the index 0 corresponds to minima, the 
index 1 corresponds to 1-saddles, the index 2 corresponds to 
2-saddles, and the index 3 corresponds to maxima. An integral 
line is a maximal curve in the domain, whose tangent aligns 
with the gradient of / at every point. The function / increases 
along the integral line and its limit points are the critical points 
of /. 

The set of all integral lines that originate at the critical 
point p together with p is called the ascending manifold of 
p. Similarly, the set of all integral lines that terminate at 
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Fig. 2: Topological cancellation of a pair of critical points 
in a 2D Morse-Smale complex, (a) Morse-Smale complex 
of the function shown in Figure 1 restricted to a 2D slice. 
Maxima are denoted by ©, saddles by ©, and minima by 
©. A pair (pi+i,qi) of critical points connected by a single 
arc is scheduled to be canceled, (b) D(pi+ 1 ) is the set of 
surviving index i critical points connected to pt+\ and A(qi) 
is the set of surviving index i + 1 critical points connected 
to qt. (c) Combinatorial realization: connect all critical points 
D(pi+ 1 ) to those in A(qi). Geometric realization: merge the 
descending manifold of pi+\ with those of critical points 
in A(qi). Merge the ascending manifold of qi with those of 
critical points in D(pi+\). 


the critical point p together with p is called the descending 
manifold of p. The ascending manifolds (similarly, the de¬ 
scending manifolds) of all critical points partition the domain. 
The ascending manifold of a critical point with index i has 
dimension n — i, where n is the dimension of the domain. Thus, 
the ascending manifold of a minimum is a three dimensional 
cell, the ascending manifold of a 1-saddle is a two dimen¬ 
sional sheet, the ascending manifold of a 2-saddle is a one 
dimensional arc (see Figure 1), and the ascending manifold 
of a maximum is equal to the maximum. The converse is 
true for the descending manifold i.e., the descending manifold 
of a critical point with index i has dimension i. The Morse- 
Smale complex is a partition of the domain into cells formed 
by the collection of integral lines that share a common source 
and a common destination. The function / is called a Morse- 
Smale function if the ascending and descending manifolds of 
all pairs of critical points intersect only transversally i.e., if 
the ascending and descending manifolds of two critical points 
intersect, then the intersection has dimension exactly equal 
to the difference in the indices of the two critical points. 
The critical points, referred to as nodes , together with the 1- 
manifolds that connect them, referred to as arcs , form the 
1-skeleton of the Morse-Smale complex, which is referred to 
as the combinatorial structure of the Morse-Smale complex. 


2.2 Morse-Smale complex simplification 

The Morse-Smale complex may be simplified by repeated 
application of the topological cancellation procedure. Topo¬ 
logical cancellation eliminates a pair of critical points in 
the Morse-Smale complex connected by a single arc. The 
resulting complex is representative of a smoother version of 
/, which may be obtained by a local smoothing operation 
within an infinitesimal neighborhood of the arc. The canceled 
pair of points are no longer critical in the smoother version 
of /. By repeatedly applying topological cancellations, one 
obtains a simpler Morse-Smale complex, where undesirable 
pairs of critical points are eliminated. This simplified Morse- 
Smale complex is representative of the function / with several 
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smoothing operations applied to it. Pairs of critical points that 
are multiply connected in the Morse-Smale complex, referred 
to as strangulations , cannot be simplified by topological 
cancellation. 

The cancellation procedure mandates changes to the com¬ 
binatorial structure and ascending/descending manifolds of 
critical points that survive the cancellation. These changes are 
respectively referred to as the combinatorial and geometric 
realization of the cancellation. Let pi+\ and qi denote the 
pair of critical points to be eliminated with i +1 and i being 
the respective Morse indices. The combinatorial realization 
proceeds by first removing pi + \ and qi as well as all arcs 
incident upon them. Next, new arcs are introduced between 
every surviving index i critical point that was connected 
to pi +1 and every surviving index i + 1 critical point that 
was connected to q t . The geometric realization merges the 
descending manifold of pi+\ with the descending manifold 
of each surviving index i + 1 critical point connected to qi. 
Analogously, the ascending manifold of qi is merged with the 
ascending manifold of each surviving index i critical point 
connected to Pi+\. Figure 2 illustrates both realizations of a 
topological cancellation procedure applied to the Morse-Smale 
complex of a two-dimensional slice of the function shown in 
Figure 1. 

The ordering of cancellation pairs plays a crucial role in de¬ 
termining the resulting structure of the Morse-Smale complex 
and its geometry. To simplify the Morse-Smale complex, pairs 
of singularly connected critical points having the least absolute 
difference in function value are iteratively canceled. This 
approach is equivalent to the notion of topological persistence 
[48], [51] for 2D Morse-Smale complexes, but not necessarily 
for 3D Morse-Smale complexes [57]. 

2.3 The hierarchical Morse-Smale complex 

A sequence of cancellations results in a hierarchical se¬ 
quence of Morse-Smale complexes MSCq,MSC\ . ....MSC n , 
where each Morse-Smale complex is a simpler version of 
the preceding Morse-Smale complex containing fewer critical 
points. Morse-Smale complex MSQ is said to be coarser than 
MSCj if i > j and finer if i < j. The version index i enumerates 
the Morse-Smale complexes in the hierarchy. Each non-zero 
version of the Morse-Smale complex, MSCi , is associated with 
the absolute difference in function value, q, of the pair of 
critical points canceled in the preceding version, MSCi- 1- As 
each iteration selects the pair of critical points with the least 
absolute difference in function value, the sequence of q’s is 
monotonically increasing i.e. (to = 0) <t\ < t n . For 

completeness of the sequence, to is set to zero. Figure 3 
illustrates a hierarchy of Morse-Smale complexes of a 2D 
equivalent of the function shown in Figure 1. 

It is not necessary to explicitly store all versions of the 
hierarchy of Morse-Smale complexes. Instead, the combi¬ 
natorial representation of only MSCo is computed initially. 
Subsequently, MSCi can be obtained from a finer MSQ-\ by 
performing topological cancellation (see Section 2.2). Analo¬ 
gously, MSCi is obtained from a coarser MSQ+ \ by applying 
the inverse operation of a topological cancellation. 



Fig. 3: Hierarchical Morse-Smale complex. A family of 
Morse-Smale complexes generated by iteratively canceling 
pairs of critical points. MSCo is the Morse-Smale complex 
of a 2D equivalent of the function shown in Figure 1. It is 
simplified to generate a coarser version, MSC\ , by canceling 
a pair of critical points (cyan) connected by a single arc 
and having least absolute difference in function value t \. 
Successive versions MSCi are computed similarly by selecting 
arcs so that to = 0 <t\ < ... < t$. 

3 Methodology 

Exploring the filamentary patterns of the cosmic web is 
challenging because of the large range of the spatial scales 
and density range it exhibits. A proper characterization should 
also account for the hierarchical nature of structures, which 
adds considerable challenges to the task. Though there exist 
different notions of filaments, the primary evidence relied upon 
for extraction and analysis is most often visual. It is therefore 
not surprising that structure finding methods often visually 
verify results by superimposing the extracted structures upon 
visualizations of the density field or the particle distribution. 
However the visualization plays a role only after structure ex¬ 
traction process in these methods. We differ in this respect by 
providing the capability to interact with the structure finding 
procedure and extract structures that are visually relevant. To 
accomplish such a visual exploration framework, a succinct 
model of filament definition, an efficient representation of 
hierarchical structures, and an appropriate query mechanism 
that supports the extraction of these structures are paramount. 
The following exposition details our framework on these 
terms. 

3.1 Density estimation and filament modeling 

Cosmological simulations are N-body particle experiments that 
simulate structure formation and evolution by tracing positions 
of the particles under the influence of physical laws. In the 
observational reality, the information about structures in the 
cosmos comes through observing the galaxies. The galaxies 
can be treated as particles also for the purpose of analysis in 
the context of large scale structures. 

The input to Felix is the logarithm of the density scalar 
field on the given 3D domain of interest. The domain could 
be 3D structured grids or tetrahedral meshes, with the density 
specified on the vertices of the grid/mesh. We find that the 
logarithm of the density field, instead of the density field 
itself, resolves the structures with more visual clarity. This 
has also been independently established in an earlier work 
[33]. Additionally, the input may be specified as a distribution 
of particles within a 3D region of interest. This could be 
a snapshot from a cosmological simulation, or galaxies in 
real observational data. We use the Delaunay tessellation field 
estimator (DTFE) [58], [59] to estimate the density of the input 
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Fig. 4: (a) A scatter plot of the function values of the canceled 
critical point pairs for the function shown in Figure 1. A 
2-saddle-maximum pair is the only pair that is far removed 
from the diagonal. This corresponds to cancellation of the 2- 
saddle with a maximum that represents one of the Gaussians 
in Figure 1. Other pairs close to the diagonal represent 
insignificant features that manifest due to the added Gaussian 
noise as well as sampling noise, (b) A scatter plot of the 
function values of the canceled critical point pairs for the 
Voronoi-Kinematic dataset B (see Section 4.1). No discernible 
separation of points is seen, though there are many points 
that are far removed from the diagonal. Thus, no clear global 
simplification threshold may be used for filament extraction. 

particles in the 3D region of interest. This procedure begins 
by computing the periodic Delaunay triangulation[60] on the 
points (simulation particles or galaxies). Next, the density at 
each vertex of the triangulation is estimated by the inverse 
of the volume of the tetrahedra incident upon it. Finally, 
the density is linearly interpolated onto the edges, faces, and 
tetrahedra of the Delaunay triangulation to yield a piecewise 
linear density function on the domain. The periodic Delaunay 
triangulation, computed by the DTFE procedure, is used to 
represent the domain. 

The Morse-Smale complex of the logarithm of the density 
field is computed. Filaments are modeled as the ascending 
manifolds of 2-saddles of the Morse-Smale complex. These 
arcs represent paths of steepest descent from the two maxima 
merging at the 2-saddle. This 2-saddle represents the lowest 
density point along the arcs connecting the two maxima. A 
schematic illustration of this is presented in Figure 1. There 
are many algorithms available in the literature to compute the 
3D Morse-Smale complex. The algorithms are primarily based 
on either the quasi Morse-Smale complex formulation [61], 
[62] or Forman’s [63] discrete Morse theory[64], [65], [66], 
[67]. We use a parallel algorithm based on the latter approach 
[67] resulting in fast computation even for large datasets. 

The density field is rarely smooth and several local maxima 
obscure a view of the larger scale behavior of the density 
field. This is especially true if the density field is computed 
on the raw particle distribution, where the density field tends 
to be spiky and with a lot of fluctuations in the high density 
cluster-like regions. The Morse-Smale complex is simplified 
by iteratively canceling pairs of singularly connected critical 
points with least absolute difference in function value to 
generate a hierarchy of Morse-Smale complexes. 

In most applications, a specific version of the Morse-Smale 


complex from the hierarchy is chosen based on a perceptibly 
clear separation of noise and features. One way to choose such 
a threshold separating noise and feature is by using a scatter 
plot of the function values of canceled critical point pairs 
(see Figure 4) where the lower function value among the pair 
corresponds to the x-coordinate and the higher function value 
corresponds to the y-coordinate. In datasets where topological 
features are well separated (see Figure 4a), pairs representing 
significant features appear far away and isolated from the 
diagonal. In such cases, the coarsest Morse-Smale complex 
version wherein the insignificant pairs are removed is selected 
for feature analysis/extraction. However, this strategy is not 
easily applicable to cosmology datasets (see Figure 4b). A well 
defined separation is rarely discernible, though there are many 
scatter points that are far removed from the diagonal. Hence, 
we drop the assumption that we must work with a specific 
version of the Morse-Smale complex. Instead, we query for 
features across all Morse-Smale complexes in the hierarchy, 
as discussed in the following sub-section. 

3.2 Density range based filament selection 

Cosmic filaments exhibit a large range of variation in their 
density characteristics. Indeed, one expects filaments to be 
present both in void like regions and between cluster like 
regions. While strong dense filaments in between clusters 
define the spine of the Cosmic Web, in the hierarchically 
evolving mass distribution we encounter a wide spectrum of 
ever more tenuous filaments on smaller mass scales. Small 
filaments define the directions of mass inflow into galaxies, 
and form a crucial component in the formation of galaxies 
[68]. Even more tenuous are the systems of filaments stretch¬ 
ing over the hollows of voids, often conspicuously aligned 
along the direction defined by neighbouring superstructures. 
The understanding of this network is tightly related to the 
issue of the “missing” dwarf galaxies in voids [69]. While this 
illustrates the complexity of the multiscale filigree of filaments 
in the Cosmic Web, we follow a strategy in which we focus 
our attention on specific aspects and details of the cosmic 
web. Dependent on the identity of objects and structures of 
interest, we wish to be able to zoom in on to the corresponding 
filamentary network. This is largely dependent on the mass 
scales of the objects involved, and the density values of the 
corresponding filament generating density peaks [70], [37]. 

Following this rationale, we translate this strategy into the 
use of queries that depend on the density properties of interest. 
Specifically, we query for filaments by specifying the density 
range [Mb,M e \ of the clusters they connect (the maxima at 
the end points), as well as the density range [Sb,S e \ of the 
lowest point along the connecting path (the density range of 
the 2-saddles). Figure 5 conceptually illustrates the character¬ 
ization of filaments using density ranges, where density along 
filaments varies significantly necessitating simplification. 

Algorithm Select2Saddles lists the algorithm to process 
such a query. The algorithm accepts, together with the com¬ 
binatorial Morse-Smale complex MSC , the density ranges of 
2-saddles [Sb,S e \ and maxima [Mb , M e \ as input (the subscripts 
b and e denote the beginning and ending of each density range 
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Fig. 5: (a) Filaments are modeled as the ascending paths of 
2-saddles connecting two extrema. The 2-saddles are filtered 
based on the range constraints [M b ,M e \ and [S b ,S e ] on the 
highest and lowest values respectively along the ascending 
paths. The highest values along the 2-saddle’s ascending 
manifold are at extrema and the lowest value is at the 2-saddle. 
The function along the paths needs to be simplified as it is 
rarely smooth. In the illustration, a simplification threshold of 
t reveals a filament with appropriate density characteristics. 
However, imposing such a threshold uniformly will cause 
another filament (b) having the required density characteristics 
to be destroyed. It is therefore necessary to extract filaments by 
querying all Morse-Smale complexes within a given hierarchy. 


respectively). The algorithm returns a list of 2-saddles that 
satisfy the above criteria together with the maximal Morse- 
Smale complex version in which they do so. 

Algorithm Select2S addles ([S b ,S e ], [M b ,M e \ ) 

1: Sver = Empty Map 

2: S = {s\s is a 2-saddle, S b < f(s) < S e } 

3: for all s e S do 

4: v a = Max{i\s is not simplified in MSQ} 

5: v b = Max{i\ti <M\j— f(s)} 

6: v c = Max{i\s connects distinct Maxima in MSQ} 

7: Sver[s] = Min(y a ,Vb,v c ) 

8: Sort S by Sver 
9: Ssel = Empty Set 
10: for all s e S do 
11: i = Sver[s] 

12: m a ,mb = Maxima connected to s in MSQ 

13: if Mb < f (m a ), f (mb) <M e then 

14: insert (s,Sver[s]) in Ssel 

15: return Ssel 


The algorithm begins by creating a list S of 2-saddles that 
have their function value in the given 2-saddle range [S b ,S e \. 
Then, for each 2-saddle in S , a Morse-Smale complex version 
in which it possibly connects two maxima within [M b ,M e \ is 
computed. The appropriate version is given by the minimum 
of three version indices v a ,v b and v c . 

The version index v a is the finest Morse-Smale complex 
version in which the 2-saddle s survives. In other words, s is 
canceled in MSQ a+ \ but not in MSQ a . This is pre-computed 
by examining the cancellation sequence. The version index 
Vb corresponds to the last Morse-Smale complex version at 
which the 2-saddle s connects two maxima, both with function 
value less than M e . This is possible because in successive 


versions of the Morse-Smale complex, the maxima connected 
to a 2-saddle via the same arc form an increasing sequence in 
terms of their function value. Thus, in the version where the 
absolute difference in function value of the last canceled pair 
is less than M b — f(s ), the 2-saddle s still possibly connects 
two maxima with function value less than M e . The version 
index v c is the last Morse-Smale complex version at which 
the 2-saddle s separates distinct maxima. In other words, it 
is not a strangulation in MSQ c . As a consequence of the 
cancellation preconditions, once a strangulation is created by a 
2-saddle, it may be destroyed only by canceling the 2-saddle 
with a 1-saddle. Thus there exists a maximal version index 
v c after which the 2-saddle remains connected to a single 
maximum. The version index v c is — 1 when the 2-saddle is a 
strangulation in the initial Morse-Smale complex. In this case, 
the 2-saddle is not considered in further steps and is removed 
from S. Again, this is easily pre-computed for each 2-saddle 
by examining the cancellation sequence. 

The 2-saddles in the set S are sorted based on their version 
indices. This is done to optimize switching between the 
required Morse-Smale complex versions. Next, each 2-saddle 
s in S is checked to see if it separates two maxima within 
the maxima density range [Mb , M e \. The list of 2-saddles that 
fulfill all of the above criteria is returned together with the 
associated version number of each 2-saddle. This above list of 
2-saddles is used to extract the filament geometry. Specifically, 
the ascending manifold of each 2-saddle is extracted from the 
corresponding version of the Morse-Smale complex. This may 
be done efficiently using the cancellation merge DAG data 
structure discussed by Gyulassy et al. [71]. 

In some situations, it is desirable to perform some simpli¬ 
fication to eliminate Poisson noise introduced due to meshing 
the domain. In these cases, a global simplification specif¬ 
ically for noise elimination, can be optionally introduced. 
Specifically, Algorithm Select2Saddles returns only those 
2-saddles that survive in hierarchical Morse-Smale complex 
versions above a specified threshold T s , where T s is specified 
as a normalized fraction of the range of log-density values 
(normalized to [0,1]). T s is set to 0.0 unless specifically men¬ 
tioned. Similar to the inputs of Algorithm Select2Saddles, 
T s may be updated during run-time. 

3.3 Volume Rendering 

We use volume visualizations of the density field to aid 
selection of parameters for Algorithm Select2Saddles. 
The geometry of the selected filaments using Algorithm Se- 
lect2Saddles is superimposed upon a volume rendering of 
the density field. Based on the visualization of the extracted 
filaments and the density volume rendering, the parameters 
may be adjusted so that the structures correspond with the 
density volume rendering. Figure 6 shows an example of 
the overlay of the volume visualization with the selected 
structures. Direct volume rendering of the density is often 
not effective for visualization because of the formation of 
clusters at multiple scales. Furthermore, these clusters are 
often spatially far removed from the features of interest. In 
the supplemental document, we describe an application of 
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our framework to suppress the opacity of regions spatially far 
removed from filamentary features of interest, which leads to 
a feature based volume rendering enhancement. 

4 Model Description 

In this section, we briefly describe the models used to test our 
filament detection routine. These are the Voronoi evolution 
models and ACDM cosmological simulations. 

The Voronoi models provide us with vital quantitative in¬ 
formation on the sensitivity of Felix to anisotropic filamentary 
patterns in the galaxy distribution. To this end it is of key 
importance that the Voronoi models have an a-priori known 
population fraction in different morphological elements: clus¬ 
ters, filaments, walls and voids. This makes them perfect test 
models for evaluating success and failure rates of the various 
identification methods. 

Although they involve filaments with a broad distribution 
of densities, the Voronoi models do not incorporate the multi¬ 
scale web-like patterns seen in realistic cosmological scenar¬ 
ios. To assess this aspect of the cosmic mass distribution, 
we turn to simulations of structure formation in the stan¬ 
dard ACDM cosmology. Implicitly these include all relevant 
physical and dynamical processes of the evolving cosmic 
dark matter distribution. However, as we have no control 
over all aspects of the emerging mass distribution in ACDM 
simulations, for testing purposes they are not as informative 
as Voronoi models 4 

4.1 Voronoi evolution models 

The Voronoi evolution models are a class of heuristic models 
for cellular distributions of galaxies that mimics the evolution 
of the Megaparsec universe towards a weblike pattern. They 
use Voronoi tessellations as a template for distribution of 
matter and related galaxy population [76], [77], [74], and its 
subsequent evolution. 

In these models, one begins by fixing an underlying Voronoi 
skeleton, defined by a small set of randomly distributed 
nuclei in the simulation box. One then superposes a set of N 
randomly distributed particles on this skeleton. The resulting 
spatial distribution of particles in the model is obtained by 
projecting the initially random distribution of particles on to 
the faces, edges, and nodes of the Voronoi tessellation. This 
results in a pattern in which one can distinguish four structural 
components: field particles located in the interior of Voronoi 
cells, wall particles within and around the Voronoi faces, 
filament particles within and around the Voronoi edges and 
cluster particles within and around the Voronoi nodes. 

One particular class of Voronoi clustering models are the 
Voronoi kinematic models, which seek to approximate the 
dynamical evolution of the large scale cosmic mass distri¬ 
bution. These models involve a continuous flow of galaxies 
towards the nearest wall, along a filament at the wall’s edge, 
and subsequently towards the final destination, a vertex of the 

4. In the literature, several studies use simplistic models using Voronoi 
tessellations. The models we use here are considerably more sophisticated, 
and represent a rather realistic depiction of the cosmic web in void-dominated 
cosmologies, see e.g. [72], [73], [74], [75]. 



cell 

wall 

filament 

cluster 

A 

29.88% 

43.57% 

22.20% 

4.33% 

B 

9.82% 

32.13% 

38.62% 

19.42% 

C 

3.5% 

16.50% 

28.70% 

51.30% 


TABLE 1: The relative abundance of particles in each struc¬ 
tural element throughout the course of evolution. 

Voronoi tessellation. This motion is regulated by the increase 
of mean distance between the galaxies, an expression of void 
expansion and evacuation as a function of time. 

The Voronoi models used in our experiments have 262,144 
particles distributed along the vertices, edges, faces and cells 
of the Voronoi skeleton in a box of side-length 200/z -1 Mpc. 
The skeleton is generated by 32 randomly placed nuclei in the 
box. For the least evolved stage, most of the particles are in the 
cells, while for the most evolved stage, most of the particles are 
located in and around clusters. Table 1 presents the percentage 
distribution of particles in the various structural elements, as 
it changes with time. Stage 1 (dataset A) is the least evolved, 
while Stage 3 (dataset C) is the most evolved. The particles 
in and around the nodes, edges, and walls are Gaussian 
distributed around these elements, characterized by a thickness 
scale Rf which is the standard deviation of the distribution. For 
our models, Rf = 2/i _1 Mpc. In the supplemental document, we 
present the volume renderings of the density distribution. 

4.2 ACDM cosmological simulations 

The ACDM simulations are fully physical models that trace 
the distribution and evolution of dark matter in the universe 
based on current understanding of real physical laws. Dark 
matter is the gravitationally dominant matter component in the 
Universe and constitutes the major fraction of matter. As it is 
known to only interact gravitationally, modeling the behavior 
of dark matter is computationally fast and efficient. Such dark 
matter simulations form one of the principal tools towards 
understanding the evolution of the matter distribution in the 
Universe. 

The cosmological simulations that we used follow the 
standard ACDM cosmology. In this model, the matter content 
of this Universe is dominated by collisionless cold dark 
matter (CDM) particles. The biggest contribution to the energy 
content of this Universe comes from dark energy, in the form 
of the cosmological constant A (see [78]), which drives its 
accelerated expansion at the current epoch. 

To present the results of our visual exploration framework, 
we use the Cosmogrid simulations [79]. It is a suite of 
simulations in a box of size 21 h~ l Mpc, each differing in the 
number of particles. The particular simulation we use for our 
study comprises of 512 3 particles. This is a relatively small 
scale in the context of the cosmic web. The mass resolution 
achieved is 8.21 x 10 6 M0. The initial conditions are setup 
at z = 65 using the Zel’dovich approximation [80]. The log- 
density field is available on a 128 x 128 x 128 structured grid. 

Particularly characteristic in the evolving mass distribution 
of the Cosmogrid simulation is the large central under-density, 
surrounded by a range of smaller voids near its outer edge. 
In combination with its extremely high spatial resolution and 
state-of-the-art dynamic range, this renders the Cosmogrid 
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simulation uniquely suited as a testbed for a case study of the 
internal structure of voids. It was precisely this circumstance 
that formed the rationale behind its exploitation in a previous 
study of the formation of dark halos along tenuous void 
filaments [81]. 

5 Results and Discussion 

In this section, we demonstrate and discuss the salient features 
and potential applications of Felix. First, we evaluate the fila¬ 
ments extracted using Felix and compare with those extracted 
using MMF, SpineWeb, and DisPerSE using the Voronoi 
kinematic datasets. Next, we present a visual exploration of the 
filaments in the Cosmogrid simulation. The methods described 
above were implemented and tested on a computer system 
with an Intel Xeon(R) 2.0 GHz CPU and 8GB of RAM. In 
the supplemental document, we briefly discuss the structure 
finders against which we compare Felix. 

5.1 Filaments in the Voronoi model: a comparison 

Here we present an analysis of the filamentary structures 
extracted using Felix, and compare the results with the tech¬ 
niques detailed above. The comparison study concerns the 
analysis results obtained for the set of heuristic Voronoi 
evolution models described above. Since they are input pa¬ 
rameters, in these models the classification of galaxies as 
void, wall, filament, and cluster are known a-priori. Following 
the application of one of the detection techniques we may 
then examine the validity and authenticity of the extracted 
structures by direct comparison with the true identity of a 
galaxy. 

For the comparison study, we define two measures. One 
quantifies the true detection rate of a method, the other the 
false identifications. We classify all galaxies within a distance 
d from the extracted structures to be filament and cluster 
particles and the others to be void and wall galaxies. For a 
given set of structures extracted from a given dataset and a 
distance d , the true positive classification rate TPj is defined 
as 

Tp # filament and cluster galaxies correctly classified 

# filament and cluster galaxies 

Similarly, the false positive classification rate FP^ is defined 
as 

pp # filament and cluster galaxies incorrectly classified 

# filament and cluster galaxies 

A large separation between these two measures indicates 
good discriminatory power of the classifier, and thus the 
proximity of relevant galaxies to the extracted structures. 

As we discuss in more detail below, the Felix’s true and false 
detection rates are comparable, and in some situations better, 
than those obtained by DisPerSE, SpineWeb and Nexus/MMF. 
A brief description of SpineWeb, and Nexus/MMF is provided 
in the supplemental document. As the DisPerSE methodology 
is closely related to Felix, in the following paragraphs, we 
briefly describe it and contrast it against Felix. 


Felix: Figure 6 shows the extracted filaments for the 
Voronoi kinematic datasets A, B , and C using Felix. The input 
density range parameters for Algorithm Select2Saddles 
are selected interactively, using the visualization information 
of the procedure. Each update is accomplished within 2- 
3 seconds. This enables interactive visual feedback so that 
parameters may be adjusted further in subsequent iterative 
steps. 

Comparison A: Felix and DisPerSE: DisPerSE [52], 
[53] is a closely related structure finder that also uses the 
Morse-Smale complex of the logarithm of the density field. 
It simplifies the Morse-Smale complex using Topological 
Persistence. The function values are normalized by the rms 
of the density field fluctuation with respect to the mean, and 
the significance level for simplification is quoted in this unit. 
Felix is closely related to DisPerSE as both use the Morse- 
Smale complex of the log-density field and involve feature 
extraction from it. DisPerSE defines significant features as 
only those that remain unsimplified using the user defined 
significance threshold. It ignores the density range charac¬ 
teristics of the extracted features. A significant consequence 
is that filaments within void-like regions and cluster like 
regions are ignored/simplified away. If they are retained, then 
the mixing of features causes visual clutter. Furthermore, the 
significance parameter selection is a fixed constant and visual 
interaction plays no role in its selection. In contrast, given the 
ubiquity of filaments in various density regimes, Felix allows 
for density ranged based probes into filaments, within clusters 
and voids. Furthermore, the visual interactive aspect allows 
for user engagement in parameter selection, which is crucial 
for the set of features identified. Another difference is that 
Felix uses simplification only for noise removal and not feature 
identification. 

In this experiment, we demonstrate the consequences of not 
correlating the density characteristics for filament extraction. 
Specifically, we demonstrate that the filaments extracted using 
Felix are more spatially proximal to filament and cluster 
particles in the Voronoi Kinematic simulation. Furthermore, 
we show that tuning the significance parameter is not a 
sufficient mechanism to extract the desired filaments in this 
dataset. In the next experiment, we demonstrate the exploration 
of filaments within high-density cluster like regions and low- 
density void like regions. Such a delineation, coupled with the 
visual exploration process, is not possible using DisPerSE. 

The recovery and failure rates for the Voronoi kinematic 
models A, B , and C are shown in Figure 8. The TPj and 
FPj of the Felix filaments are plotted in the top-left panel, 
the ones for DisPerSE can be found in the top right-hand 
panel (Id significance threshold), bottom left-hand panel (3d 
significance threshold) and bottom-right panel (5d significance 
threshold). 

Felix shows good recovery rates for all datasets, particularly 
around d = 3/z -1 Mpc. For the least evolved configuration A, 
and for locations where structures are least distinct, Felix 
still obtains moderately good recovery rates. The results for 
DisPerSE with simplification thresholds Id and 3d are com¬ 
parable: at short distances the true detection rate is slightly 
lower than that of Felix, while at larger distances it performs 
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(a) (b) (c) 

Fig. 6: Filamentary structures extracted from datasets A, B, and C of the Voronoi evolution time-series using Felix with 
parameters for Algorithm Select2Saddles set as follows: (a) [Sb,S e \ = [10 -L6 ,oo] and [Mb,M e \ = [10 _01 ,oo], (b,c) [Sb,S e \ = 
[10- L6 ,oo] and [Mb,M e \ = [10°,oo]. Filaments are shown as orange tubes along with a volume rendering of the log-density 
field. The dense knot like structures show filaments within cluster-like regions. 



Fig. 7: Filamentary structures extracted from datasets A, B and C of the Voronoi evolution time-series using DisPerSE with 
significance level of 5 a. Filaments are shown as orange tubes along with a volume rendering of the log-density field. The 
inset pictures show identified filaments that are within wall-like and void-like regions of the Voronoi kinematic datasets. 



Fig. 8: Classification recovery rates for Voronoi kinematic 
datasets A, B, and C. Top left: Felix; top right: DisPerSE, 
1 a significance level; bottom left: DisPerSE, 3 a significance 
level, and (bottom right) DisPerSE with a 5cr significance 
level. False positive rates greater than 1.0 are clipped and 
respective values are shown. 


marginally better. The situation is slightly different in the case 
of DisPerSE with a 5(7 simplification threshold. In the case 
of the more strongly evolved B and C datasets, DisPerSE and 
Felix have similar false detection rates FPj, while the true 
detection rates TPd of Felix are consistently higher. 

For the lesser evolved datasets A and B , the false detection 
rates FPj for DisPerSE quickly increase to rather large values. 
For simplification thresholds of la and 3 a it even surpasses 
values of unity. This may indicate that in certain circumstances 
an automatic detection of filaments from the Morse-Smale 
complex runs the risk of over-determining the population of 
filaments, even after considerable simplification. While the 
problem is not so acute in the most evolved stage C, where 
the morphologies are well separated, direct simplification 
strategies may not always succeed in properly classifying all 
filament, wall and void regions in the more moderately evolved 
stages A and B. 

Figure 7 shows the filaments in the three Voronoi models 
detected by DisPerSE, with a simplification threshold of 5a. 
In comparison with the structures in Figure 6, the knot like 
structures present in clusters are absent. This leads to the 
cluster particles being far away from the filament end points, 
and thus the reduced TPj rates of DisPerSE. In contrast, 
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Dataset 


TP d 

FP d 

A 

Felix(d = 3 h : Mpc) 

0.67 

0.36 

A 

DisPerSE (5 o,d = 3/z -1 Mpc) 

0.59 

0.31 

B 

Felix(d = 3/i _1 Mpc) 

0.84 

0.13 

B 

DlsPerSE{ 5cr,rf = S/r'Mpc) 

0.69 

0.11 


Nexus/MMF 

0.85 

0.13 

C 

Felix(d = 3/z _1 Mpc) 

0.90 

0.05 

C 

DisPerSE (5 o,d = 3/z -1 Mpc) 

0.78 

0.05 


SpineWeb 

0.87 

0.10 


TABLE 2: Recovery rates of galaxies within 3/z _1 Mpc of 
structures extracted using Felix compared with Nexus/MMF, 
SpineWeb, and DisPerSE (5cr significance level). 

Felix’s ranged query allows us to retain only the filaments 
in cluster like regions and those that connect these cluster 
like regions, leading to better TP a rates. Also, many filaments 
found by DisPerSE are within the wall-like and void-like 
regions of the Voronoi Kinematic simulation. Some examples 
are highlighted using insets in Figure 7. Again these are 
filtered out by Felix’s ranged query, which is not directly 
possible in DisPerSE. The inclusion of such structures in 
DisPerSE leads to its higher FPj rates. 

These findings suggest that a structure identification strategy 
based on a direct simplification procedure of the Morse- 
Smale complex should be applied with care to the density 
regimes being studied. Specifically, the superior classification 
rate profiles confirm that, using Felix, we can easily extract 
filaments that are spatially more proximal to the cluster and 
filament particles in the Voronoi Kinematic datasets. This 
issue is not easily addressed by the significance threshold of 
DisPerSE. In contrast, Felix provides an intutive density based 
handle to extract the desired features. Also, the possibility 
of having FPj value larger than 1 in extreme situations is 
indicative of over-detection of filaments. This is potentially 
cumbersome for the analysis of genuine cosmological simu¬ 
lations and observational surveys. In more complex realistic 
circumstances, cosmic structure involves features over a wide 
range of densities and scales and structural morphologies that 
are not as well separated as in the simpler Voronoi models. 

Comparison B: Felix and SpineWeb: For the com¬ 
parison of the Felix and Spineweb [70], we use a Voronoi 
Evolution model realization that is comparable to the advanced 
state of dataset C. We use the test result reported in [70] 
with respect to the model that has a similar percentage of 
particles in the four morphological features. In these model 
realizations, the clusters, filaments and walls have a Gaussian 
density profile with a scale of R g = l/z -1 Mpc. The spine has 
an effective width d = 2 R g as the identified structures are 
thickened by 1 voxel, with a size of R g = l/z -1 Mpc. 

For this configuration, Aragon-Calvo et al. [70] report true 
and false detection rates of TPd = 0.87 and FPd = 0.10 (see 
Table 2). Felix attains the same recovery rate of TPd = 0.87 at 
a smaller distance d = 3/* _1 Mpc = 1.5 R g . For the same config¬ 
uration, the failure rate parameter, FPd = 0.05, is comparable 
to that reported for Spineweb. By comparison, at d = 2 R g , for 
Felix the recovery rates are TPd = 0.93 and FPd = 0.07. 

In summary, these results appear to suggest that Felix 
performs as well as SpineWeb. 


Comparison C: Felix and Nexus/MMF: Since 
Nexus/MMF concerns a sophisticated formalism based on a 
scale space analysis, the parameters of detection do not have 
a direct correspondence with topology based techniques like 
DisPerSE, SpineWeb, and Felix. 

We use dataset B for a comparison with Nexus/MMF. This 
dataset is similar to the least evolved dataset used by [31] 
in an evaluation of the MMF, the original density field based 
Nexus/MMF implementation. For similar values of the detec¬ 
tion rate TPd, both Felix and MMF have identical failure rates 
FPj. This indicates that both procedures have a comparable 
detection behavior. 

5.2 Filament Exploration 

In this section we discuss the application of Felix to explore 
different classes of filaments from cosmological simulations. 
The ability to filter filaments on a combination of morpho¬ 
logical and density properties is helpful in situations where 
we wish to focus on, for example, the properties of galaxies 
residing in filaments in low-density void regions or in the high- 
density outskirts of clusters. This ability of Felix to identify 
a specified population of intravoid filaments or cluster inflow 
channels provides us with a microscopic instrument that allows 
a detailed and systematic exploration of the fine structure in 
the hierarchy of cosmic structure. Felix is able to zoom in on 
such regions and delineate their detailed infrastructure. 

As the criteria for the identification of filaments and other 
web-like features still differ substantially between the various 
available techniques, the visual interaction aspect of Felix 
is a major practical asset in obtaining a proper user-defined 
selection of filaments. To this end, we may also point out 
that available automatic detection techniques may produce 
significant spurious results, which may substantially influence 
the results of targeted studies as the one illustrated here. A 
telling example of this has been discussed in the previous 
section. 

Figures 9a and 9c present volume renderings of a 3D region 
of the Cosmogrid dataset ranging from z-coordinates 69 to 
105. The bounding box of the dataset is 128 x 128 x 128. 
This region is selected as it contains a large void like region 
surrounded by a large number of high density regions. The 
transfer function opacities have been adjusted to highlight 
the filament like structures in cluster-like and void-like re¬ 
gions respectively. The hierarchical Morse-Smale complex 
computation and the filament selection are executed on the 
entire dataset. The resulting filaments and volume rendering 
are clipped to the above mentioned region of interest. The 
filament selection parameters of algorithm Select2Saddles 
are adjusted interactively, via a visualization step (see ac¬ 
companying video). Each selection query takes approximately 
1 second to process. The subsequent extraction of filament 
geometry depends on the number of selected 2-saddles. This 
takes approximately 4 seconds. Thus, the query framework 
may be used to interactively change parameters and visually 
correlate the set of extracted features with the underlying 
density distribution. 

The filament selections obtained following the application 
of the interactive procedure are illustrated Figures 9b and 
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Fig. 9: Exploring filaments in high-density cluster-like environments and low-density void like regions in the Cosmogrid dataset. 
The selected region, comprising of a 3D region of z-slices from 69 through 105, contains filamentary structures in both types of 
environments. Volume rendering of the density field with opacity adjusted to highlight filaments in (a) high density cluster-like 
regions, and (c) low density void-like regions, (b) Filaments within high density cluster-like regions extracted with parameters 
[Sb,S e ] = [10°, 10 9 * 6 ] and [M^ , M e \ = [10 23 ,10 9 6 ]. (d) Filaments within low density void-like regions extracted with parameters 
[‘S'&jSej = [10 -2 - 5 ,10° 5 ] and = [10 a5 ,10 3 * 5 ]. For both sets of filaments, the value of T s is set to 0.05. 


9d. Each shows filaments in a different environmental density 
regime. Figure 9b shows filaments that exist near and within 
the high density cluster like regions. These are the filaments 
that form the spine of the cosmic web. Figure 9d shows 
filaments within void-like regions. The combination of density 
criteria and interactive visualization enables us to zoom in on 
this system of intra-void filaments. They are the faint residuals 
of the smaller-scale filaments that constituted the spine of the 
cosmic web at earlier cosmic epochs, and as such represent 
a direct manifestation of the hierarchical buildup of cosmic 
structure. At the current epoch, the intra-void filaments appear 
to define a different pattern than the prominent filamentary 
bridges between clusters of galaxies. As a result of the tidal 
influence of surrounding large-scale mass concentrations they 
are conspicuously aligned along a direction correlated with the 
main axis of the embedding void. 

6 Conclusion 

We have presented a Topology based Framework, named 
Felix, to probe filament structure in the large scale universe. 
The framework is particularly designed to probe filamentary 


structures in different density regimes, and optimally preserve 
structural detail in regimes of interest. While other cosmic 
structure analysis tools do not include a facility to select web¬ 
like features according to tailor-made aspects and characteris¬ 
tics, this is precisely the mission of the Felix procedure. We 
directed Felix towards a case study of the the filamentary in¬ 
frastructure and architecture of cosmic voids and demonstrated 
that it successfully extracts the network of tenuous filaments 
pervading their interior [72], [73], [74], [75] . 

In an accompanying study, we plan to exploit the Felix 
facility to study the physical characteristics of the extracted 
samples of intra-void filaments. This also involves their halo 
and subhalo population, their gas content, and the relation of 
these with the embedding voids and surrounding large-scale 
mass distribution. This will be of key importance towards 
understanding the formation and evolution of void galaxies 
[82], [83], [84] and specifically that of the issue of the missing 
dwarf galaxies [69]. In addition, following the recognition that 
void architecture represents a potentially sensitive probe of 
dark energy and dark matter and a keen test of modified gravity 
theories [9], [10], [11], [12], [13], [14], [15], [16], the filament 
samples extracted by Felix will be subjected to a systematic 
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study of their dependence on cosmological parameters. 

As an immediate extension, we plan to use Felix with other 
scalar fields such as the tidal force field. Another possible 
direction is the visualization and analysis of the hierarchy of 
voids, walls, and filaments in cosmological datasets. Interac¬ 
tive visual exploration of these intricate structural networks 
remains a challenging and largely unexplored problem of 
major significance. 
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